#############################################################################################
########### calculate exposure inex for each species ########################################
#############################################################################################

#load packages:
library(R.utils)
library(raster)
library(sp)
library(sf)
library(stringr)
library(circular)
library(plotrix)
library(grDevices)

#Main working dir:
MyDir<-paste("/home/jc217070/Maxent",sep="")  #define working directory on HPC

#MyDir<-paste("C:/Users/AnnaF/Dropbox/temp/WHOSnakeData/")
#PlotDir<-MyDir
OccDir<-file.path(MyDir,"Maxent_SWDs/spp_swds")
MUOutDir<-paste(MyDir, "/Outputs2_Crossval/permutation.importance", sep="")
IndSppOutDir<-paste(MyDir, "/Outputs3_Final/Final_Individual_Spp", sep="")
PlotDir<-paste(MyDir,"/Maxent_Data/Plots/Plots_WHOSnakes_2020_09_01",sep="")

#set project name:
projectName<-"WHOSnakes"

#define 'for' variables:
crit<-c("permutation.importance")   #c("permutation.importance","contribution","Test.gain.without","Test.gain.with.only","AUC.without","AUC.with.only")
timecats<-c("Current","2050_Median","2050_Q90","2050_Q10","2090_Median","2090_Q90","2090_Q10")
#timecats<-c("Current","2050_Median","2090_Median")
#timecats<-c("Current")

#read species reference list:
DatSum<-read.table(file=file.path(paste(MyDir, "/Maxent_LUTables/",projectName,".csv",sep="")), header=TRUE, sep=",") #load data summary for all species
#DatSum<-read.table(file=file.path(paste(MyDir, "/2_Maxent_LUTables/",projectName,".csv",sep="")), header=TRUE, sep=",") #load data summary for all species
DatSum<-DatSum[DatSum$Listing_comment!="not listed",]
spp<-DatSum$gen_sp_subtax[DatSum$Listing_comment!="not listed"]
sppA<-DatSum$Accepted_Species[DatSum$Listing_comment!="not listed"]
#spSTART<-which (spp=="Micrurus_mertensi")

DatSumMod<-read.table(file=paste(MyDir, "/Maxent_LUTables/",projectName,"_SupMatFinFin.csv",sep=""), header=TRUE, sep=",") #load data summary for all species
#
#DatSumMod<-read.table(file=paste(MyDir, "/2_Maxent_LUTables/",projectName,"_SupMatFinFin.csv",sep=""), header=TRUE, sep=",") #load data summary for all species
# DatSumMod$RangesizeSQKM_Current[DatSumMod$RangesizeSQKM_Current==0]<-0.01
# DatSumMod$SuitabilitySum_Current[DatSumMod$SuitabilitySum_Current==0]<-0.01
# DatSumMod$Exposure_Current[DatSumMod$Exposure_Current==0]<-0.01
# DatSumMod$Impact_Current[DatSumMod$Impact_Current==0]<-0.01






bg<-raster(paste(PlotDir,"/Cop_FAPAR_mean_2010-2019_WW1km.asc",sep=""), crs="+init=epsg:4326")
myext<-extent(-170,179,-58,72)
bg<-crop(bg, myext)
GLOBraster<-bg
values(GLOBraster)<-ifelse(values(GLOBraster)!=-Inf,0,-Inf) #raster showing all land as 0, ocean as NA

GLOBrasterB<-GLOBraster
values(GLOBrasterB)<-ifelse(values(GLOBrasterB)==0,1,values(GLOBrasterB)) #raster showing all land as 1, ocean as NA
GLOBrasterC<-GLOBraster
GLOBrasterC[GLOBrasterC == 0] <- NA






worldmap  <- st_read(paste(PlotDir,"/worldmap/WHO_countries.shp",sep=""))
WHOpolies <- st_read(paste(MyDir,"/WHOSnakesShapes/SSN_CTY_CAT_BRG_FEA_IN_Mar2023.shp",sep=""))
SDMregions   <- st_read(paste(PlotDir,"/worldmap/SDMRegions_WGS84.shp",sep=""))

#make region data frame columns
########################################################################################################################################################################################################################
#basic info:
myregions <- unique(SDMregions$Short_Name)     #vector of region names
myregionsLong <- unique(SDMregions$Region_Nam)     #vector of region names
# region.tab <- data.frame(myregions,myregionsLong) #empty vector to add data for regions into
# names(region.tab)<-c("region_name_short","region_name_long")
# region.tab$n_spp                       <- NA
# ########################################################################################################################################################################################################################
# #columns for range size changes:
# #2050
# region.tab$n_RangesizeSQKM_contraction5to50P_2050_Median<-NA 
# region.tab$n_RangesizeSQKM_expansion5to50P_2050_Median <-NA
# region.tab$n_RangesizeSQKM_contraction50plusP_2050_Median <-NA
# region.tab$n_RangesizeSQKM_expansion50plusP_2050_Median<-NA
# region.tab$n_RangesizeSQKM_stable_2050_Median<-NA
# #2090
# region.tab$n_RangesizeSQKM_contraction5to50P_2090_Median<-NA
# region.tab$n_RangesizeSQKM_expansion5to50P_2090_Median <-NA
# region.tab$n_RangesizeSQKM_contraction50plusP_2090_Median <-NA
# region.tab$n_RangesizeSQKM_expansion50plusP_2090_Median <-NA
# region.tab$n_RangesizeSQKM_stable_2090_Median <-NA
# ########################################################################################################################################################################################################################
# #columns for suitability changes:
# #2050
# region.tab$n_SuitabilitySum_contraction5to50P_2050_Median<-NA
# region.tab$n_SuitabilitySum_expansion5to50P_2050_Median<-NA
# region.tab$n_SuitabilitySum_contraction50plusP_2050_Median<-NA
# region.tab$n_SuitabilitySum_expansion50plusP_2050_Median <-NA
# region.tab$n_SuitabilitySum_stable_2050_Median <-NA
# #2090
# region.tab$n_SuitabilitySum_contraction5to50P_2090_Median<-NA
# region.tab$n_SuitabilitySum_expansion5to50P_2090_Median<-NA
# region.tab$n_SuitabilitySum_contraction50plusP_2090_Median <-NA
# region.tab$n_SuitabilitySum_expansion50plusP_2090_Median<-NA
# region.tab$n_SuitabilitySum_stable_2090_Median<-NA
# ########################################################################################################################################################################################################################
# #columns for exposure changes:
# #2050
# region.tab$n_Exposure_contraction5to50P_2050_Median<-NA
# region.tab$n_Exposure_expansion5to50P_2050_Median<-NA
# region.tab$n_Exposure_contraction50plusP_2050_Median<-NA
# region.tab$n_Exposure_expansion50plusP_2050_Median<-NA
# region.tab$n_Exposure_stable_2050_Median<-NA
# #2090
# region.tab$n_Exposure_contraction5to50P_2090_Median <-NA 
# region.tab$n_Exposure_expansion5to50P_2090_Median<-NA
# region.tab$n_Exposure_contraction50plusP_2090_Median <-NA
# region.tab$n_Exposure_expansion50plusP_2090_Median<-NA
# region.tab$n_Exposure_stable_2090_Median<-NA
# ########################################################################################################################################################################################################################
# #columns for shift vectors:
# #2050
# region.tab$med_ShiftDist_2050_Median<-NA
# region.tab$med_ShiftDir_2050_Median<-NA
# region.tab$mean_ShiftDist_2050_Median<-NA
# region.tab$mean_ShiftDir_2050_Median<-NA
# region.tab$q10_ShiftDist_2050_Median<-NA
# region.tab$q10_ShiftDir_2050_Median <-NA
# region.tab$q90_ShiftDist_2050_Median<-NA
# region.tab$q90_ShiftDir_2050_Median <-NA
# region.tab$med_ShiftDist_2050_Q10 <-NA
# region.tab$med_ShiftDir_2050_Q10 <-NA
# region.tab$med_ShiftDist_2050_Q90 <-NA
# region.tab$med_ShiftDir_2050_Q90<-NA
# #2090
# region.tab$med_ShiftDist_2090_Median<-NA
# region.tab$med_ShiftDir_2090_Median<-NA
# region.tab$mean_ShiftDist_2090_Median<-NA
# region.tab$mean_ShiftDir_2090_Median<-NA
# region.tab$q10_ShiftDist_2090_Median<-NA
# region.tab$q10_ShiftDir_2090_Median <-NA
# region.tab$q90_ShiftDist_2090_Median<-NA
# region.tab$q90_ShiftDir_2090_Median <-NA
# region.tab$med_ShiftDist_2090_Q10 <-NA
# region.tab$med_ShiftDir_2090_Q10<-NA 
# region.tab$med_ShiftDist_2090_Q90 <-NA 
# region.tab$med_ShiftDir_2090_Q90<-NA
# ########################################################################################################################################################################################################################
# 
# 
# 
# 
# for (region in myregions){
#  
#   print(paste("filling out data table for region",region))
#   n<- which(region.tab$region_name_short==region)
#   regionpoly <- subset(SDMregions, SDMregions$Short_Name == region)     #polygon for this region
#   sf_use_s2(FALSE)
#   mypoliesTF<-st_intersects(WHOpolies,regionpoly,sparse=FALSE) #species polygons overlapping with this region polygon
#   WHOpolies$TF<-mypoliesTF
#   mypolies<-subset(WHOpolies, WHOpolies$TF=="TRUE")
#   myspp<- sort(gsub("_"," ",unique(mypolies$Scientific)))    #vector of species present in this region
#   myspp<- myspp[myspp %in% gsub("_"," ",DatSum$gen_sp_subtax)] #safety, don't use any polygons for species that are not listed
#   region.tab$n_spp[n] <- length(myspp)
#   
#   ########################################################################################################################################################################################################################
#   #columns for range size changes:
#   #2050
#   region.tab$n_RangesizeSQKM_contraction5to50P_2050_Median[n]    <- length(DatSumMod$RangesizeSQKM_CHANGE_2050_Median[DatSumMod$RangesizeSQKM_CHANGE_2050_Median/DatSumMod$RangesizeSQKM_Current   <  (-0.05) &
#                                                                                                                         DatSumMod$RangesizeSQKM_CHANGE_2050_Median/DatSumMod$RangesizeSQKM_Current >  (-0.5)  &
#                                                                                                                         DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.tab$n_RangesizeSQKM_expansion5to50P_2050_Median[n]      <- length(DatSumMod$RangesizeSQKM_CHANGE_2050_Median[DatSumMod$RangesizeSQKM_CHANGE_2050_Median/DatSumMod$RangesizeSQKM_Current   >    0.05  &
#                                                                                                                         DatSumMod$RangesizeSQKM_CHANGE_2050_Median/DatSumMod$RangesizeSQKM_Current <    0.5   &
#                                                                                                                         DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.tab$n_RangesizeSQKM_contraction50plusP_2050_Median[n]   <- length(DatSumMod$RangesizeSQKM_CHANGE_2050_Median[DatSumMod$RangesizeSQKM_CHANGE_2050_Median/DatSumMod$RangesizeSQKM_Current   <= (-0.5)  &
#                                                                                                                         DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.tab$n_RangesizeSQKM_expansion50plusP_2050_Median[n]     <- length(DatSumMod$RangesizeSQKM_CHANGE_2050_Median[DatSumMod$RangesizeSQKM_CHANGE_2050_Median/DatSumMod$RangesizeSQKM_Current   >=   0.5   &
#                                                                                                                         DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.tab$n_RangesizeSQKM_stable_2050_Median[n]               <- length(DatSumMod$RangesizeSQKM_CHANGE_2050_Median[DatSumMod$RangesizeSQKM_CHANGE_2050_Median/DatSumMod$RangesizeSQKM_Current   <=    0.05  &
#                                                                                                                         DatSumMod$RangesizeSQKM_CHANGE_2050_Median/DatSumMod$RangesizeSQKM_Current >=  (-0.05) &
#                                                                                                                         DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   #2090
#   region.tab$n_RangesizeSQKM_contraction5to50P_2090_Median[n]    <- length(DatSumMod$RangesizeSQKM_CHANGE_2090_Median[DatSumMod$RangesizeSQKM_CHANGE_2090_Median/DatSumMod$RangesizeSQKM_Current   <  (-0.05) &
#                                                                                                                         DatSumMod$RangesizeSQKM_CHANGE_2090_Median/DatSumMod$RangesizeSQKM_Current >  (-0.5)  &
#                                                                                                                         DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.tab$n_RangesizeSQKM_expansion5to50P_2090_Median[n]      <- length(DatSumMod$RangesizeSQKM_CHANGE_2090_Median[DatSumMod$RangesizeSQKM_CHANGE_2090_Median/DatSumMod$RangesizeSQKM_Current   >    0.05  &
#                                                                                                                         DatSumMod$RangesizeSQKM_CHANGE_2090_Median/DatSumMod$RangesizeSQKM_Current <    0.5   &
#                                                                                                                         DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.tab$n_RangesizeSQKM_contraction50plusP_2090_Median[n]   <- length(DatSumMod$RangesizeSQKM_CHANGE_2090_Median[DatSumMod$RangesizeSQKM_CHANGE_2090_Median/DatSumMod$RangesizeSQKM_Current   <= (-0.5)  &
#                                                                                                                         DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.tab$n_RangesizeSQKM_expansion50plusP_2090_Median[n]     <- length(DatSumMod$RangesizeSQKM_CHANGE_2090_Median[DatSumMod$RangesizeSQKM_CHANGE_2090_Median/DatSumMod$RangesizeSQKM_Current   >=   0.5   &
#                                                                                                                         DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.tab$n_RangesizeSQKM_stable_2090_Median[n]               <- length(DatSumMod$RangesizeSQKM_CHANGE_2090_Median[DatSumMod$RangesizeSQKM_CHANGE_2090_Median/DatSumMod$RangesizeSQKM_Current   <=    0.05  &
#                                                                                                                         DatSumMod$RangesizeSQKM_CHANGE_2090_Median/DatSumMod$RangesizeSQKM_Current >=  (-0.05) &
#                                                                                                                         DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   ########################################################################################################################################################################################################################
#   #columns for suitability changes:
#   #2050
#   region.tab$n_SuitabilitySum_contraction5to50P_2050_Median[n]    <- length(DatSumMod$SuitabilitySum_CHANGE_2050_Median[DatSumMod$SuitabilitySum_CHANGE_2050_Median/DatSumMod$SuitabilitySum_Current   <  (-0.05) &
#                                                                                                                         DatSumMod$SuitabilitySum_CHANGE_2050_Median/DatSumMod$SuitabilitySum_Current   >  (-0.5)  &
#                                                                                                                         DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.tab$n_SuitabilitySum_expansion5to50P_2050_Median[n]      <- length(DatSumMod$SuitabilitySum_CHANGE_2050_Median[DatSumMod$SuitabilitySum_CHANGE_2050_Median/DatSumMod$SuitabilitySum_Current   >    0.05  &
#                                                                                                                         DatSumMod$SuitabilitySum_CHANGE_2050_Median/DatSumMod$SuitabilitySum_Current   <    0.5   &
#                                                                                                                         DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.tab$n_SuitabilitySum_contraction50plusP_2050_Median[n]   <- length(DatSumMod$SuitabilitySum_CHANGE_2050_Median[DatSumMod$SuitabilitySum_CHANGE_2050_Median/DatSumMod$SuitabilitySum_Current   <= (-0.5)  &
#                                                                                                                         DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.tab$n_SuitabilitySum_expansion50plusP_2050_Median[n]     <- length(DatSumMod$SuitabilitySum_CHANGE_2050_Median[DatSumMod$SuitabilitySum_CHANGE_2050_Median/DatSumMod$SuitabilitySum_Current   >=   0.5   &
#                                                                                                                         DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.tab$n_SuitabilitySum_stable_2050_Median[n]               <- length(DatSumMod$SuitabilitySum_CHANGE_2050_Median[DatSumMod$SuitabilitySum_CHANGE_2050_Median/DatSumMod$SuitabilitySum_Current   <=    0.05  &
#                                                                                                                         DatSumMod$SuitabilitySum_CHANGE_2050_Median/DatSumMod$SuitabilitySum_Current   >=  (-0.05) &
#                                                                                                                         DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   #2090
#   region.tab$n_SuitabilitySum_contraction5to50P_2090_Median[n]    <- length(DatSumMod$SuitabilitySum_CHANGE_2090_Median[DatSumMod$SuitabilitySum_CHANGE_2090_Median/DatSumMod$SuitabilitySum_Current   <  (-0.05) &
#                                                                                                                         DatSumMod$SuitabilitySum_CHANGE_2090_Median/DatSumMod$SuitabilitySum_Current   >  (-0.5)  &
#                                                                                                                         DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.tab$n_SuitabilitySum_expansion5to50P_2090_Median[n]      <- length(DatSumMod$SuitabilitySum_CHANGE_2090_Median[DatSumMod$SuitabilitySum_CHANGE_2090_Median/DatSumMod$SuitabilitySum_Current   >    0.05  &
#                                                                                                                         DatSumMod$SuitabilitySum_CHANGE_2090_Median/DatSumMod$SuitabilitySum_Current   <    0.5   &
#                                                                                                                         DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.tab$n_SuitabilitySum_contraction50plusP_2090_Median[n]   <- length(DatSumMod$SuitabilitySum_CHANGE_2090_Median[DatSumMod$SuitabilitySum_CHANGE_2090_Median/DatSumMod$SuitabilitySum_Current   <= (-0.5)  &
#                                                                                                                         DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.tab$n_SuitabilitySum_expansion50plusP_2090_Median[n]     <- length(DatSumMod$SuitabilitySum_CHANGE_2090_Median[DatSumMod$SuitabilitySum_CHANGE_2090_Median/DatSumMod$SuitabilitySum_Current   >=   0.5   &
#                                                                                                                         DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.tab$n_SuitabilitySum_stable_2090_Median[n]               <- length(DatSumMod$SuitabilitySum_CHANGE_2090_Median[DatSumMod$SuitabilitySum_CHANGE_2090_Median/DatSumMod$SuitabilitySum_Current   <=    0.05  &
#                                                                                                                         DatSumMod$SuitabilitySum_CHANGE_2090_Median/DatSumMod$SuitabilitySum_Current   >=  (-0.05) &
#                                                                                                                         DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   ########################################################################################################################################################################################################################
#   #columns for exposure changes:
#   #2050
#   region.tab$n_Exposure_contraction5to50P_2050_Median[n]    <- length(DatSumMod$Exposure_CHANGE_2050_Median[DatSumMod$Exposure_CHANGE_2050_Median/DatSumMod$Exposure_Current                        <  (-0.05) &
#                                                                                                                           DatSumMod$Exposure_CHANGE_2050_Median/DatSumMod$Exposure_Current          >  (-0.5)  &
#                                                                                                                           DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.tab$n_Exposure_expansion5to50P_2050_Median[n]      <- length(DatSumMod$Exposure_CHANGE_2050_Median[DatSumMod$Exposure_CHANGE_2050_Median/DatSumMod$Exposure_Current                        >    0.05  &
#                                                                                                                           DatSumMod$Exposure_CHANGE_2050_Median/DatSumMod$Exposure_Current          <    0.5   &
#                                                                                                                           DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.tab$n_Exposure_contraction50plusP_2050_Median[n]   <- length(DatSumMod$Exposure_CHANGE_2050_Median[DatSumMod$Exposure_CHANGE_2050_Median/DatSumMod$Exposure_Current                        <= (-0.5)  &
#                                                                                                                           DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.tab$n_Exposure_expansion50plusP_2050_Median[n]     <- length(DatSumMod$Exposure_CHANGE_2050_Median[DatSumMod$Exposure_CHANGE_2050_Median/DatSumMod$Exposure_Current                        >=   0.5   &
#                                                                                                                           DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.tab$n_Exposure_stable_2050_Median[n]               <- length(DatSumMod$Exposure_CHANGE_2050_Median[DatSumMod$Exposure_CHANGE_2050_Median/DatSumMod$Exposure_Current                        <=    0.05  &
#                                                                                                                           DatSumMod$Exposure_CHANGE_2050_Median/DatSumMod$Exposure_Current          >=  (-0.05) &
#                                                                                                                           DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   #2090
#   region.tab$n_Exposure_contraction5to50P_2090_Median[n]    <- length(DatSumMod$Exposure_CHANGE_2090_Median[DatSumMod$Exposure_CHANGE_2090_Median/DatSumMod$Exposure_Current                        <  (-0.05) &
#                                                                                                                           DatSumMod$Exposure_CHANGE_2090_Median/DatSumMod$Exposure_Current          >  (-0.5)  &
#                                                                                                                           DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.tab$n_Exposure_expansion5to50P_2090_Median[n]      <- length(DatSumMod$Exposure_CHANGE_2090_Median[DatSumMod$Exposure_CHANGE_2090_Median/DatSumMod$Exposure_Current                        >    0.05  &
#                                                                                                                           DatSumMod$Exposure_CHANGE_2090_Median/DatSumMod$Exposure_Current          <    0.5   &
#                                                                                                                           DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.tab$n_Exposure_contraction50plusP_2090_Median[n]   <- length(DatSumMod$Exposure_CHANGE_2090_Median[DatSumMod$Exposure_CHANGE_2090_Median/DatSumMod$Exposure_Current                        <= (-0.5)  &
#                                                                                                                           DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.tab$n_Exposure_expansion50plusP_2090_Median[n]     <- length(DatSumMod$Exposure_CHANGE_2090_Median[DatSumMod$Exposure_CHANGE_2090_Median/DatSumMod$Exposure_Current                        >=   0.5   &
#                                                                                                                           DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.tab$n_Exposure_stable_2090_Median[n]               <- length(DatSumMod$Exposure_CHANGE_2090_Median[DatSumMod$Exposure_CHANGE_2090_Median/DatSumMod$Exposure_Current                        <=    0.05  &
#                                                                                                                           DatSumMod$Exposure_CHANGE_2090_Median/DatSumMod$Exposure_Current          >=  (-0.05) &
#                                                                                                                           DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   ########################################################################################################################################################################################################################
#   #columns for shift vectors:
#   #2050
#   region.tab$med_ShiftDist_2050_Median[n]            <- median(DatSumMod$DistM_2050_Median[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)], na.rm=TRUE)
#   region.tab$med_ShiftDir_2050_Median[n]             <- median.circular(circular(DatSumMod$DirDeg_2050_Median[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)],
#                                                                                  units = "degrees",modulo = "asis",rotation = "clock"), na.rm=TRUE)
#   region.tab$mean_ShiftDist_2050_Median[n]            <- mean(DatSumMod$DistM_2050_Median[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)], na.rm=TRUE)
#   region.tab$mean_ShiftDir_2050_Median[n]             <- mean.circular(circular(DatSumMod$DirDeg_2050_Median[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)],
#                                                                                 units = "degrees",modulo = "asis",rotation = "clock"), na.rm=TRUE)
#   region.tab$q10_ShiftDist_2050_Median[n]            <- quantile(DatSumMod$DistM_2050_Median[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)],probs=0.1, na.rm=TRUE)
#   region.tab$q10_ShiftDir_2050_Median[n]             <- quantile.circular(circular(DatSumMod$DirDeg_2050_Median[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)],
#                                                                               units = "degrees",modulo = "asis",rotation = "clock"),probs=0.1, na.rm=TRUE)
#   region.tab$q90_ShiftDist_2050_Median[n]            <- quantile(DatSumMod$DistM_2050_Median[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)],probs=0.9, na.rm=TRUE)
#   region.tab$q90_ShiftDir_2050_Median[n]             <- quantile.circular(circular(DatSumMod$DirDeg_2050_Median[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)],
#                                                                               units = "degrees",modulo = "asis",rotation = "clock"), probs=0.9,na.rm=TRUE)
#   region.tab$med_ShiftDist_2050_Q10[n]               <- median(DatSumMod$DistM_2050_Q10[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)], na.rm=TRUE)
#   region.tab$med_ShiftDir_2050_Q10[n]                <- median.circular(circular(DatSumMod$DirDeg_2050_Q10[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)],
#                                                                                  units = "degrees",modulo = "asis",rotation = "clock"), na.rm=TRUE)
#   region.tab$med_ShiftDist_2050_Q90[n]               <- median(DatSumMod$DistM_2050_Q90[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)], na.rm=TRUE)
#   region.tab$med_ShiftDir_2050_Q90[n]                <- median.circular(circular(DatSumMod$DirDeg_2050_Q90[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)],
#                                                                                units = "degrees",modulo = "asis",rotation = "clock"), na.rm=TRUE)
#   #2090
#   region.tab$med_ShiftDist_2090_Median[n]            <- median(DatSumMod$DistM_2090_Median[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)], na.rm=TRUE)
#   region.tab$med_ShiftDir_2090_Median[n]             <- median.circular(circular(DatSumMod$DirDeg_2090_Median[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)],
#                                                                                  units = "degrees",modulo = "asis",rotation = "clock"), na.rm=TRUE)
#   region.tab$mean_ShiftDist_2090_Median[n]            <- mean(DatSumMod$DistM_2090_Median[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)], na.rm=TRUE)
#   region.tab$mean_ShiftDir_2090_Median[n]             <- mean.circular(circular(DatSumMod$DirDeg_2090_Median[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)],
#                                                                                  units = "degrees",modulo = "asis",rotation = "clock"), na.rm=TRUE)
#   region.tab$q10_ShiftDist_2090_Median[n]            <- quantile(DatSumMod$DistM_2090_Median[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)],probs=0.1, na.rm=TRUE)
#   region.tab$q10_ShiftDir_2090_Median[n]             <- quantile.circular(circular(DatSumMod$DirDeg_2090_Median[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)],
#                                                                                    units = "degrees",modulo = "asis",rotation = "clock"), probs=0.1,na.rm=TRUE)
#   region.tab$q90_ShiftDist_2090_Median[n]            <- quantile(DatSumMod$DistM_2090_Median[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)],probs=0.9, na.rm=TRUE)
#   region.tab$q90_ShiftDir_2090_Median[n]             <- quantile.circular(circular(DatSumMod$DirDeg_2090_Median[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)],
#                                                                                    units = "degrees",modulo = "asis",rotation = "clock"), probs=0.9,na.rm=TRUE)
#   region.tab$med_ShiftDist_2090_Q10[n]               <- median(DatSumMod$DistM_2090_Q10[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)], na.rm=TRUE)
#   region.tab$med_ShiftDir_2090_Q10[n]                <- median.circular(circular(DatSumMod$DirDeg_2090_Q10[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)],
#                                                                                  units = "degrees",modulo = "asis",rotation = "clock"), na.rm=TRUE)
#   region.tab$med_ShiftDist_2090_Q90[n]               <- median(DatSumMod$DistM_2090_Q90[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)], na.rm=TRUE)
#   region.tab$med_ShiftDir_2090_Q90[n]                <- median.circular(circular(DatSumMod$DirDeg_2090_Q90[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)],
#                                                                                  units = "degrees",modulo = "asis",rotation = "clock"), na.rm=TRUE)
#   ########################################################################################################################################################################################################################
# }
# 
# write.table(x=region.tab, file=paste(MyDir, "/Maxent_LUTables/",projectName,"_RegionData.csv",sep=""),
#             na = "NA", row.names = FALSE, col.names = TRUE, sep=",")
# 
# 
# 
# ######################
# #plot bar graphs for regions:
# 

region.tab<-read.table(file=paste(MyDir, "/Maxent_LUTables/",projectName,"_RegionData.csv",sep=""), 
                       header=TRUE, sep=",")


#rearrange order of regions for plot:
myregions<-c(myregions[3],myregions[4],myregions[2],myregions[10],myregions[6],
              myregions[7],myregions[1],myregions[8],myregions[5],myregions[9])

png(filename=paste(PlotDir,"/Region_CCBars.png",sep=""), units="in", width=16, height=23, res=500) #png smaller file size, A3 for better resolution
par(mfrow = c(5,2), mar = c(7,6,5,1), oma=c(0,0.2,0.2,0.2), mgp=c(2.7,0.8,0))
# #mar=c(bottom, left, top, right); oma= outer margins; mpg = axis label locations c(3, 1, 0)


for (region in myregions){

  n<- which(region.tab$region_name_short==region)

  # Create matrix
  set.seed(112)
  data <- matrix(nrow=5,ncol=6)
  colnames(data) <- c("Range Size 2050","Range Size 2090","Suitability 2050","Suitability 2090","Exposure 2050","Exposure 2090")
  rownames(data) <- c("<-50%","<-5%","+/-5%",">+5%",">+50%")

  # Get data for this region:
  data[,"Range Size 2050"]<-c(region.tab$n_RangesizeSQKM_contraction50plusP_2050_Median[n],
                              region.tab$n_RangesizeSQKM_contraction5to50P_2050_Median[n],
                              region.tab$n_RangesizeSQKM_stable_2050_Median[n],
                              region.tab$n_RangesizeSQKM_expansion5to50P_2050_Median[n],
                              region.tab$n_RangesizeSQKM_expansion50plusP_2050_Median[n])/region.tab$n_spp[n]
  data[,"Range Size 2090"]<-c(region.tab$n_RangesizeSQKM_contraction50plusP_2090_Median[n],
                              region.tab$n_RangesizeSQKM_contraction5to50P_2090_Median[n],
                              region.tab$n_RangesizeSQKM_stable_2090_Median[n],
                              region.tab$n_RangesizeSQKM_expansion5to50P_2090_Median[n],
                              region.tab$n_RangesizeSQKM_expansion50plusP_2090_Median[n])/region.tab$n_spp[n]
  data[,"Suitability 2050"]<-c(region.tab$n_SuitabilitySum_contraction50plusP_2050_Median[n],
                              region.tab$n_SuitabilitySum_contraction5to50P_2050_Median[n],
                              region.tab$n_SuitabilitySum_stable_2050_Median[n],
                              region.tab$n_SuitabilitySum_expansion5to50P_2050_Median[n],
                              region.tab$n_SuitabilitySum_expansion50plusP_2050_Median[n])/region.tab$n_spp[n]
  data[,"Suitability 2090"]<-c(region.tab$n_SuitabilitySum_contraction50plusP_2090_Median[n],
                              region.tab$n_SuitabilitySum_contraction5to50P_2090_Median[n],
                              region.tab$n_SuitabilitySum_stable_2090_Median[n],
                              region.tab$n_SuitabilitySum_expansion5to50P_2090_Median[n],
                              region.tab$n_SuitabilitySum_expansion50plusP_2090_Median[n])/region.tab$n_spp[n]
  data[,"Exposure 2050"]<-c(region.tab$n_Exposure_contraction50plusP_2050_Median[n],
                               region.tab$n_Exposure_contraction5to50P_2050_Median[n],
                               region.tab$n_Exposure_stable_2050_Median[n],
                               region.tab$n_Exposure_expansion5to50P_2050_Median[n],
                               region.tab$n_Exposure_expansion50plusP_2050_Median[n])/region.tab$n_spp[n]
  data[,"Exposure 2090"]<-c(region.tab$n_Exposure_contraction50plusP_2090_Median[n],
                               region.tab$n_Exposure_contraction5to50P_2090_Median[n],
                               region.tab$n_Exposure_stable_2090_Median[n],
                               region.tab$n_Exposure_expansion5to50P_2090_Median[n],
                               region.tab$n_Exposure_expansion50plusP_2090_Median[n])/region.tab$n_spp[n]


  # Get the stacked barplot
  barplot(data,
          col=c("royalblue","lightblue1","grey90","pink","red"),
          border="white",
          space=c(0.1,0.1,0.4,0.1,0.4,0.1),
          ylim=c(0,1),
          font.axis=1, cex.axis=2,cex.names=2,
          xlab=NULL,ylab="fraction of species",cex.lab=2,
          names.arg=c("2050","2090","2050","2090","2050","2090"),
          legend.tex=FALSE,
          main=region.tab$region_name_long[n],cex.main=2)
  #args.legend = list(x = "topright",inset = c(-0.2, 0),cex=0.8)
  lines (c(-0.2,7.3),c(0.5,0.5), col="darkgrey", lwd= 2, lty=1)
  lines (c(-0.2,7.3),c(0,0), col="black", lwd= 2, lty=1)
  
  #add text for range size, suitability, etc below x axis:
  # text(0.5,0.5, labels=paste("Range Size"),pos=4, col="black", cex=0.5)
  # text(3,-1.5, labels=paste("Suitability"),pos=4, col="black", cex=0.5)
  # text(5.5,-1.5, labels=paste("Exposure"),pos=4, col="black", cex=0.5)

  mtext(text=paste("Range Size"),  font= 1, side = 1, line = 4, cex= 1.5, adj=0.08)
  mtext(text=paste("Suitability"), font= 1, side = 1, line = 4, cex= 1.5, adj=0.51)
  mtext(text=paste("SHOI"),        font= 1, side = 1, line = 4, cex= 1.5, adj=0.86)

  # text(0.05, data[1,1]/2, labels=paste(">50%\ndecrease"),pos=4, col="black", cex=0.9)
  # text(0.05, data[1,1]+data[2,1]/2, labels=paste("5% to 50%\ndecrease"),pos=4, col="black", cex=0.9)
  # text(0.05, data[1,1]+data[2,1]+data[3,1]/2, labels=paste("stable\n(0% to 5%\nchange)"),pos=4, col="black", cex=0.9)
  # text(0.05, data[1,1]+data[2,1]+data[3,1]+data[4,1]/2, labels=paste("5% to 50%\nincrease"),pos=4, col="black", cex=0.9)
  # text(0.05, data[1,1]+data[2,1]+data[3,1]+data[4,1]+data[5,1]/2, labels=paste(">50%\nincrease"),pos=4, col="black", cex=0.9)


}

dev.off()

# #############################
# #plot circular graphs for regions
# 
# 
# 
# #############################
# #plot map with species vectors, and region vectors
# 
mycols1<-c("grey90","grey85","grey80","grey75","grey70",
           "grey65","grey60","grey55","grey50","grey45",
           "grey40","grey35","grey30")
xs <- c(24,55,46,133,79,-58,114,-101,102,-99)
ys <- c(8,40,58,-25,22,-12,0,21,39,45)


png(filename=paste(PlotDir,"/regionshiftmap.png",sep=""), units="in", width=20, height=16, res=500) #png smaller file size, A3 for better resolution
par(mfrow = c(2,1), mar = c(0.5,2,0,0), oma=c(0,0,1,0), mgp=c(1.2,0.5,0))
#mar=c(bottom, left, top, right); oma= outer margins; mpg = axis label locations c(3, 1, 0)

#Plot 1: 2050
#plot background
plot(bg, col=mycols1, box=FALSE, bty="n", cex.axis=0.8, legend=FALSE, axes=FALSE,main=NULL)
plot(worldmap[5], border="grey20", lwd= 1, add=TRUE,col=NA,reset=FALSE)
plot(SDMregions[2], border="black", lwd= 3, add=TRUE,col=NA,reset=FALSE)

#plot species vectors
for (sp in spp){
  l<-which (DatSumMod$Accepted_Species==gsub("_"," ",spp))
  arrows(x0=DatSumMod$long_Current[l],
         y0=DatSumMod$lat_Current[l],
         x1=DatSumMod$long_2050_Median[l],
         y1=DatSumMod$lat_2050_Median[l],
         length= 0.08, angle=20,col="black",lty=1,lwd=1.5)
}

for (region in myregions){

  n<- which(region.tab$region_name_short==region)
  s<- which(myregions==region)

  DirDeg<-    c(region.tab$med_ShiftDir_2050_Q10[n],
                region.tab$med_ShiftDir_2050_Q90[n],
                region.tab$med_ShiftDir_2050_Median[n])
  DirDist<-   c(region.tab$med_ShiftDist_2050_Q10[n],
                region.tab$med_ShiftDist_2050_Q90[n],
                region.tab$med_ShiftDist_2050_Median[n])/1000/100*20 # meters div by 1000=km, div by 100=dd
  DirDegCirc<- circular(DirDeg, type = "directions",
                        units = "degrees",
                        modulo = "asis",
                        rotation = "clock")
  #Draw arrows with package circular:
  arrows.circular(DirDegCirc,DirDist, x0 = rep(xs[s],3), y0 = rep(ys[s],3), na.rm = FALSE,zero = (pi/2),
                  col=c("royalblue","orangered","black"),
                  lwd=c(4,4,4),lty=c(1,1,1))
}



#Plot 2: 2090
#plot background
plot(bg, col=mycols1,box=FALSE, bty="n", cex.axis=0.8, legend=FALSE, axes=FALSE,main=NULL)
plot(worldmap[5], border="grey20", lwd= 0.5, add=TRUE,col=NA,reset=FALSE)
plot(SDMregions[2], border="black", lwd= 2, add=TRUE,col=NA,reset=FALSE)

#plot species vectors
for (sp in spp){
  l<-which (DatSumMod$Accepted_Species==gsub("_"," ",spp))
  arrows(x0=DatSumMod$long_Current[l],
         y0=DatSumMod$lat_Current[l],
         x1=DatSumMod$long_2090_Median[l],
         y1=DatSumMod$lat_2090_Median[l],
         length= 0.08, angle=20,col="black",lty=1,lwd=1.5)
}

for (region in myregions){

  n<- which(region.tab$region_name_short==region)
  s<- which(myregions==region)

  DirDeg<-    c(region.tab$med_ShiftDir_2090_Q10[n],
                region.tab$med_ShiftDir_2090_Q90[n],
                region.tab$med_ShiftDir_2090_Median[n])
  DirDist<-   c(region.tab$med_ShiftDist_2090_Q10[n],
                region.tab$med_ShiftDist_2090_Q90[n],
                region.tab$med_ShiftDist_2090_Median[n])/1000/100*20 # meters div by 1000=km, div by 100=dd
  DirDegCirc<- circular(DirDeg, type = "directions",
                      units = "degrees",
                      modulo = "asis",
                      rotation = "clock")
  #Draw arrows with package circular:
    arrows.circular(DirDegCirc,DirDist, x0 = rep(xs[s],3), y0 = rep(ys[s],3), na.rm = FALSE,zero = (pi/2),
                    col=c("royalblue","orangered","black"),
                    lwd=c(4,4,4),lty=c(1,1,1))
}

dev.off()

##############################################################
# show main 10 SPECIES per REGION that have the highest current exposure index
# 
DatSumMod<-read.table(file=paste(MyDir, "/Maxent_LUTables/",projectName,"_SupMatFinFin.csv",sep=""), header=TRUE, sep=",") #load data summary for all species
DatSumMod<-DatSumMod[DatSumMod$Accepted_Species %in% sppA,]
region.tab<-read.table(file=paste(MyDir, "/Maxent_LUTables/",projectName,"_RegionData.csv",sep=""), header=TRUE, sep=",") #load data summary for all species
RegionExps<-read.table(file=paste(MyDir, "/Maxent_LUTables/PerBioRegionExposures.csv",sep=""), header=TRUE, sep=",") #load data summary for all species

png(filename=paste(PlotDir,"/Region_Main10Spp.png",sep=""), units="in", width=16, height=23, res=500) #png smaller file size, A3 for better resolution
par(mfrow = c(5,2), mar = c(13,6,6,1), oma=c(0,0.2,0.2,0.2), mgp=c(4,0.5,0))
# #mar=c(bottom, left, top, right); oma= outer margins; mpg = axis label locations c(3, 1, 0)

for (region in myregions){

  n<- which(region.tab$region_name_short==region)
  s<- which(myregions==region)

  # print(paste(Sys.time(), "extracting exposure data per species for", myregions[s], sep=" "))
  # 
  # regionpoly <- subset(SDMregions, SDMregions$Short_Name == region)     #polygon for this region
  # regionXs<-st_coordinates(regionpoly)[,1]
  # regionYs<-st_coordinates(regionpoly)[,2]
  # myextMOD<-extent(regionpoly)
  # snapraster<-crop(GLOBrasterC, myextMOD)
  # regionras  <- rasterize(x=regionpoly,y=snapraster,as.numeric(1))
  # 
  # sf_use_s2(FALSE)
  # mypoliesTF<-st_intersects(WHOpolies,regionpoly,sparse=FALSE) #species polygons overlapping with this region polygon
  # WHOpolies$TF<-mypoliesTF
  # mypolies<-subset(WHOpolies, WHOpolies$TF=="TRUE")
  # myspp<- sort(gsub("_"," ",unique(mypolies$Scientific)))    #vector of species present in this region
  # myspp<- myspp[myspp %in% gsub("_"," ",DatSum$gen_sp_subtax)] #safety, don't use any polygons for species that are not listed
  # 
  # sp.exps<-vector()
  # 
  # for (sp in myspp){
  # 
  #   print(paste("calculating exposure for",sp,"within",region))
  # 
  #   outfile<-paste(IndSppOutDir,"/",gsub(" ","_",sp),"_Current_Exp.asc",sep="")
  # 
  #   # unzip raster
  #   if(file.exists(paste(outfile,".gz", sep="")) & !file.exists(outfile)){
  #     gunzip(filename=paste(outfile,".gz", sep=""),
  #            destname=paste(outfile),
  #            remove=FALSE, overwrite=TRUE)
  #   }
  # 
  #   # read raster:
  #   Exp<-raster(paste(outfile), crs="+init=epsg:4326")
  #   #only keep raster values within region:
  #   Exp<-Exp*regionras
  # 
  #   myvals <- values(Exp)
  #   myvals<- myvals[myvals>0]
  #   myvals<- myvals[!is.na(myvals)]
  #   mysum<-sum(myvals)
  # 
  #   sp.exps<-c(sp.exps,mysum)
  # 
  #   #remove unzipped raster or zip
  #   if(file.exists(paste(outfile,".gz", sep="")) & file.exists(outfile)){
  #     file.remove(outfile)
  #   }
  #   if(!file.exists(paste(outfile,".gz", sep="")) & file.exists(outfile)){
  #     gzip(filename=paste(outfile))
  #   }
  # 
  # }
  # 
  # # DatSumModRegion<-DatSumMod[DatSumMod$Accepted_Species %in% myspp,]
  # DatSumModRegion<-data.frame(myspp,sp.exps,rep(region,length(myspp)))
  # names(DatSumModRegion)<-c("Accepted_Species","Exposure_Current","Region")
  # DatSumModRegion<-DatSumModRegion[order(DatSumModRegion$Exposure_Current,decreasing=TRUE),]
  # 
  # print(paste(Sys.time(),"exposure for",region,":"))
  # print(paste(DatSumModRegion))
  # 
  # if(s==1){
  #   RegionExps<-DatSumModRegion
  # } else{
  #   RegionExps<-rbind(RegionExps,DatSumModRegion)
  # }
  # 
  
  
  
  DatSumModRegion<-RegionExps[RegionExps$Region==region,]
  
  
  
  
  sp.exps<-DatSumModRegion$Exposure_Current[1:10]
  tot.exp<-sum(DatSumModRegion$Exposure_Current,na.rm=TRUE)
  res.exp<-tot.exp-sum(sp.exps,na.rm=TRUE)
  sp.exps<-c(sp.exps,res.exp)
  top10<-DatSumModRegion$Accepted_Species[1:10]
  top10<-gsub(" ","\n",c(top10,'other'))

  sp.exps[11]<-ifelse(sp.exps[11] > sp.exps[1],sp.exps[1],sp.exps[11])
  
  barplot(sp.exps/1000000,las=2,
          col=c(heat.colors(11,1,rev=FALSE)),
          border="grey40",
          space=c(0.1),
          ylim=c(0,max(DatSumModRegion$Exposure_Current/1000000)),
          font.axis=1,
          xlab=NULL,ylab="SHOI [millions]",
          names.arg=top10,
          legend.tex=FALSE,
          main=region.tab$region_name_long[n],
          cex.main=2,cex.axis=1.8,cex.lab=2,cex.names=1.7)
  lines (c(-0.2,12.1),c(0,0), col="black", lwd= 2, lty=1)
  
  if(sp.exps[11] >= sp.exps[1]){
    
    if(round(res.exp/1000000)>=10){
      text(x=11.2,y=(sp.exps[1]/20/1000000)*16, labels=paste(round(res.exp/1000000)),pos=4, col="black", cex=2)
    }else{
      text(x=11.3,y=(sp.exps[1]/20/1000000)*16, labels=paste(round(res.exp/1000000)),pos=4, col="black", cex=2)
    }
    
    lines (c(11,12.2),c((sp.exps[1]/20/1000000)*10,
                      (sp.exps[1]/20/1000000)*12),
           col="black", lwd= 2, lty=1)
    lines (c(11,12.2),c((sp.exps[1]/20/1000000)*12,
                      (sp.exps[1]/20/1000000)*14),
           col="black", lwd= 2, lty=1)
    
  }
  
  # write.table(x=RegionExps, file=paste(MyDir, "/Maxent_LUTables/PerBioRegionExposures_prelim.csv",sep=""),
  #             na = "NA", row.names = FALSE, col.names = TRUE, sep=",")
  # 
  # polygon(regionpoly, border = NA,
  #         col=rgb(col2rgb("black")[1],col2rgb("black")[2],col2rgb("black")[3],max=255,alpha=75))

}

dev.off()

 # write.table(x=RegionExps, file=paste(MyDir, "/Maxent_LUTables/PerBioRegionExposures.csv",sep=""),
 #             na = "NA", row.names = FALSE, col.names = TRUE, sep=",")
# 
# 



 ##############################################################
 # show ALL TAXONOMIC GROUPS per REGION that have the highest current exposure index
 # 
 DatSumMod<-read.table(file=paste(MyDir, "/Maxent_LUTables/",projectName,"_SupMatFinFin.csv",sep=""), header=TRUE, sep=",") #load data summary for all species
 DatSumMod<-DatSumMod[DatSumMod$Accepted_Species %in% sppA,]
 region.tab<-read.table(file=paste(MyDir, "/Maxent_LUTables/",projectName,"_RegionData.csv",sep=""), header=TRUE, sep=",") #load data summary for all species
 RegionExps<-read.table(file=paste(MyDir, "/Maxent_LUTables/PerBioRegionExposures.csv",sep=""), header=TRUE, sep=",") #load data summary for all species
 DatSum<-DatSum[DatSum$Accepted_Species %in% sppA,]
 mygroups<-unique(DatSum$Label)
 
 # png(filename=paste(PlotDir,"/Region_Main10SppGroups.png",sep=""), units="in", width=16, height=23, res=500) #png smaller file size, A3 for better resolution
 # par(mfrow = c(5,2), mar = c(13,6,6,1), oma=c(0,0.2,0.2,0.2), mgp=c(4,0.5,0))
 # # #mar=c(bottom, left, top, right); oma= outer margins; mpg = axis label locations c(3, 1, 0)
 
 for (region in myregions){
   
   n<- which(region.tab$region_name_short==region)
   s<- which(myregions==region)
   ThisRegionExps<-RegionExps[RegionExps$Region==region,]
   
   group.exps<-vector()
   region.groups<-vector()
   
      for (group in mygroups){
     
     m<-which(mygroups==group)
     
     myspp<-DatSum$Accepted_Species[DatSum$Label==group] #list all species in this group
     myrows<-which(ThisRegionExps$Accepted_Species %in% myspp)      #get row numbers of region exposures that belong to species from this group
     
     if(length(myrows)>=1){
       
     mysum<-sum(ThisRegionExps$Exposure_Current[myrows])
     group.exps<-c(group.exps,mysum)
     region.groups<-c(region.groups,group)
     
     }
   }
   

   DatSumModRegion<-data.frame(region.groups,group.exps,rep(region,length(region.groups)))
   names(DatSumModRegion)<-c("TaxonomicGroup","Exposure_Current","Region")
   DatSumModRegion<-DatSumModRegion[order(DatSumModRegion$Exposure_Current,decreasing=TRUE),]
   
   print(paste(Sys.time(),"exposure for",region,":"))
   print(paste(DatSumModRegion))
   
   if(s==1){
     RegionExpsGroups<-DatSumModRegion
   } else{
     RegionExpsGroups<-rbind(RegionExpsGroups,DatSumModRegion)
   }
   
   # group.exps<-DatSumModRegion$Exposure_Current[1:10]
   # tot.exp<-sum(DatSumModRegion$Exposure_Current,na.rm=TRUE)
   # res.exp<-tot.exp-sum(group.exps,na.rm=TRUE)
   # group.exps<-c(group.exps,res.exp)
   # top10<-DatSumModRegion$MergeShort[1:10]
   # top10<-gsub(" ","\n",c(top10,'other'))
   # 
   # group.exps[11]<-ifelse(group.exps[11] > group.exps[1],group.exps[1],group.exps[11])
   # 
   # barplot(group.exps/1000000,las=2,
   #         col=c(heat.colors(11,1,rev=FALSE)),
   #         border="grey40",
   #         space=c(0.1),
   #         ylim=c(0,max(DatSumModRegion$Exposure_Current/1000000)),
   #         font.axis=1,
   #         xlab=NULL,ylab="SHOI [millions]",
   #         names.arg=top10,
   #         legend.tex=FALSE,
   #         main=region.tab$region_name_long[n],
   #         cex.main=2,cex.axis=1.8,cex.lab=2,cex.names=1.7)
   # lines (c(-0.2,12.1),c(0,0), col="black", lwd= 2, lty=1)
   # 
   # if(group.exps[11] >= group.exps[1]){
   #   
   #   if(round(res.exp/1000000)>=10){
   #     text(x=11.2,y=(group.exps[1]/20/1000000)*16, labels=paste(round(res.exp/1000000)),pos=4, col="black", cex=2)
   #   }else{
   #     text(x=11.3,y=(group.exps[1]/20/1000000)*16, labels=paste(round(res.exp/1000000)),pos=4, col="black", cex=2)
   #   }
   #   
   #   lines (c(11,12.2),c((group.exps[1]/20/1000000)*10,
   #                       (group.exps[1]/20/1000000)*12),
   #          col="black", lwd= 2, lty=1)
   #   lines (c(11,12.2),c((group.exps[1]/20/1000000)*12,
   #                       (group.exps[1]/20/1000000)*14),
   #          col="black", lwd= 2, lty=1)
   #   
   # }
   # 
   
   
   group.exps<-DatSumModRegion$Exposure_Current
   top10<-DatSumModRegion$TaxonomicGroup
   top10<-gsub(" ","\n",c(top10))
   
   png(filename=paste(PlotDir,"/",region,"RankedSppGroups.png",sep=""), units="in", width=25, height=15, res=500) #png smaller file size, A3 for better resolution
   par(mfrow = c(1,1), mar = c(15,6,6,1), oma=c(0,0.2,0.2,0.2), mgp=c(4,0.5,0))
   
   barplot(group.exps/1000000,las=2,
           col=c(heat.colors(length(top10),1,rev=FALSE)),
           border="grey40",
           space=c(0.1),
           ylim=c(0,max(DatSumModRegion$Exposure_Current/1000000)),
           font.axis=1,
           xlab=NULL,ylab="SHOI [millions]",
           names.arg=top10,
           legend.tex=FALSE,
           main=region.tab$region_name_long[n],
           cex.main=2,cex.axis=1.8,cex.lab=2,cex.names=1.7)
   lines (c(-0.2,length(top10)+1.1),c(0,0), col="black", lwd= 2, lty=1)
   
   dev.off()
   
   write.table(x=RegionExpsGroups, file=paste(MyDir, "/Maxent_LUTables/PerBioRegionExposuresSppGroups_prelim.csv",sep=""),
               na = "NA", row.names = FALSE, col.names = TRUE, sep=",")
   
 }
 
 
 write.table(x=RegionExpsGroups, file=paste(MyDir, "/Maxent_LUTables/PerBioRegionExposuresSppGroups.csv",sep=""),
             na = "NA", row.names = FALSE, col.names = TRUE, sep=",")
 
 
 
 
 
 
 
 ##############################################################
 # show ALL TAXONOMIC GROUPS per COUNTRY that have the highest current exposure index
 # 
 DatSumMod<-read.table(file=paste(MyDir, "/Maxent_LUTables/",projectName,"_SupMatFinFin.csv",sep=""), header=TRUE, sep=",") #load data summary for all species
 DatSumMod<-DatSumMod[DatSumMod$Accepted_Species %in% sppA,]
 region.tab<-read.table(file=paste(MyDir, "/Maxent_LUTables/",projectName,"_RegionData.csv",sep=""), header=TRUE, sep=",") #load data summary for all species
 RegionExps<-read.table(file=paste(MyDir, "/Maxent_LUTables/PerCountryExposures.csv",sep=""), header=TRUE, sep=",") #load data summary for all species
 DatSum<-DatSum[DatSum$Accepted_Species %in% sppA,]

 #mycountries <- gsub("C<U+00F4>te D<U+2019>Ivoire","Cote d'Ivoire",worldmap$CODE_CTY)   #vector of region names
 mycountries<-unique(RegionExps$CountryCode)
 mycountrynames<-unique(RegionExps$CountryName)
 
 
 
 for (country in mycountries){
   
   
   n<- which(mycountries==country)
   #countryName<- gsub("C<U+00F4>te D<U+2019>Ivoire","Cote d'Ivoire",worldmap$NAME_CTY[n])
   s<- which(mycountries==country)
   countryname<-mycountrynames[n]
   ThisRegionExps<-RegionExps[RegionExps$CountryCode==country,]
   
   group.exps<-vector()
   region.groups<-vector()
   
   
   # print(paste(Sys.time(), "extracting exposure data per species for", countryName, sep=" "))
   # 
   # regionpoly <- subset(worldmap, worldmap$CODE_CTY == country)     #polygon for this region
   # # regionXs<-st_coordinates(regionpoly)[,1]
   # # regionYs<-st_coordinates(regionpoly)[,2]
   # myextMOD<-extent(regionpoly)
   # 
   # sf_use_s2(FALSE)
   # mypoliesTF<-st_intersects(WHOpolies,regionpoly,sparse=FALSE) #species polygons overlapping with this region polygon
   # WHOpolies$TF<-mypoliesTF
   # mypolies<-subset(WHOpolies, WHOpolies$TF=="TRUE")
   # myspp<- sort(gsub("_"," ",unique(mypolies$Scientific)))    #vector of species present in this region
   # myspp<- myspp[myspp %in% gsub("_"," ",sppA)] #safety, don't use any polygons for species that are not listed
   # 
   # if(length(myspp)>=1){
   #   
   #   snapraster<-crop(GLOBrasterC, myextMOD)
   #   regionras  <- rasterize(x=regionpoly,y=snapraster,as.numeric(1))
   #   
   #   sp.exps<-vector()
   #   
   #   for (sp in myspp){
   #     
   #     print(paste(Sys.time(),"calculating exposure for",sp,"within",country))
   #     
   #     outfile<-paste(IndSppOutDir,"/",gsub(" ","_",sp),"_Current_Exp.asc",sep="")
   #     
   #     # unzip raster
   #     if(file.exists(paste(outfile,".gz", sep="")) & !file.exists(outfile)){
   #       gunzip(filename=paste(outfile,".gz", sep=""),
   #              destname=paste(outfile),
   #              remove=FALSE, overwrite=TRUE)
   #     }
   #     
   #     # read raster:
   #     Exp<-raster(paste(outfile), crs="+init=epsg:4326")
   #     #only keep raster values within region:
   #     Exp<-Exp*regionras
   #     
   #     myvals <- values(Exp)
   #     myvals<- myvals[myvals>0]
   #     myvals<- myvals[!is.na(myvals)]
   #     mysum<-sum(myvals)
   #     
   #     sp.exps<-c(sp.exps,mysum)
   #     
   #     #remove unzipped raster or zip
   #     if(file.exists(paste(outfile,".gz", sep="")) & file.exists(outfile)){
   #       file.remove(outfile)
   #     }
   #     if(!file.exists(paste(outfile,".gz", sep="")) & file.exists(outfile)){
   #       gzip(filename=paste(outfile))
   #     }
   #     
   #   }
   #   
   #   # DatSumModRegion<-DatSumMod[DatSumMod$Accepted_Species %in% myspp,]
   #   DatSumModRegion<-data.frame(myspp,sp.exps)
   #   names(DatSumModRegion)<-c("Accepted_Species","Exposure_Current")
   #   DatSumModRegion<-DatSumModRegion[order(DatSumModRegion$Exposure_Current,decreasing=TRUE),]
   #   DatSumModRegion$CountryCode<-paste(country)
   #   DatSumModRegion$CountryName<-paste(countryName)
   #   
   #   print(paste(Sys.time(),"exposure for",country,":"))
   #   print(paste(DatSumModRegion))
   #   
   #   if(n==1){
   #     CountryExps<-DatSumModRegion
   #   } else{
   #     CountryExps<-read.table(file=paste(MyDir, "/Maxent_LUTables/PerCountryExposures.csv",sep=""), header=TRUE, sep=",") #load data summary for all species
   #     CountryExps<-rbind(CountryExps,DatSumModRegion)
   #   }
   #   
   #   write.table(x=CountryExps, file=paste(MyDir, "/Maxent_LUTables/PerCountryExposures.csv",sep=""),
   #               na = "NA", row.names = FALSE, col.names = TRUE, sep=",")
   #   
     # sp.exps<-DatSumModRegion$Exposure_Current[1:10]
     # tot.exp<-sum(DatSumModRegion$Exposure_Current,na.rm=TRUE)
     # res.exp<-tot.exp-sum(sp.exps,na.rm=TRUE)
     # sp.exps<-c(sp.exps,res.exp)
     # top10<-DatSumModRegion$Accepted_Species[1:10]
     # top10<-gsub(" ","\n",c(top10,'other'))
     # 
     # png(filename=paste(PlotDir,"/",country,"Main10Spp.png",sep=""), units="in", width=16, height=16, res=500) #png smaller file size, A3 for better resolution
     # par(mfrow = c(1,1), mar = c(10,10,5,1), oma=c(0,0.5,2,0.5), mgp=c(5,0.5,0))
     # # #mar=c(bottom, left, top, right); oma= outer margins; mpg = axis label locations c(3, 1, 0)
     # 
     # 
     # barplot(sp.exps,las=2,
     #         col=c(heat.colors(11,1,rev=FALSE)),
     #         border="white",
     #         space=c(0.1),
     #         ylim=c(0,max(DatSumModRegion$Exposure_Current)),
     #         font.axis=1, cex.axis=1.2,
     #         xlab=NULL,ylab="EI",cex.lab=1.2,
     #         names.arg=top10,
     #         legend.tex=FALSE,
     #         main=countryName,cex.main=1.8)
     # 
   
     DatSumModRegion<-ThisRegionExps
     
     print(paste(Sys.time(),"exposure for",country,":"))
     print(paste(DatSumModRegion))
     
     group.exps<-DatSumModRegion$Exposure_Current
     #tot.exp<-sum(DatSumModRegion$Exposure_Current,na.rm=TRUE)
     #res.exp<-tot.exp-sum(group.exps,na.rm=TRUE)
     #group.exps<-c(group.exps,res.exp)
     top10<-DatSumModRegion$Accepted_Species
     top10<-gsub(" ","\n",c(top10))
     
     #group.exps[11]<-ifelse(group.exps[11] > group.exps[1],group.exps[1],group.exps[11])
     
     png(filename=paste(PlotDir,"/",country,"_RankedSpp.png",sep=""), units="in", width=45, height=15, res=500) #png smaller file size, A3 for better resolution
     par(mfrow = c(1,1), mar = c(15,6,6,1), oma=c(0,0.2,0.2,0.2), mgp=c(4,0.5,0))
     
     barplot(group.exps/1000000,las=2,
             col=c(heat.colors(length(top10),1,rev=FALSE)),
             border="grey40",
             space=c(0.1),
             ylim=c(0,max(DatSumModRegion$Exposure_Current/1000000)),
             font.axis=1,
             xlab=NULL,ylab="SHOI [millions]",
             names.arg=top10,
             legend.tex=FALSE,
             main=countryname,
             cex.main=2,cex.axis=1.8,cex.lab=2,cex.names=1.7)
     lines (c(-0.2,length(top10)+1.1),c(0,0), col="black", lwd= 2, lty=1)
     
     dev.off()
     
   #   removeTmpFiles(h=0.1) #removes all temporary files older than 30min, then read hotspot map back in to continue process
   #   
   # } else{
   #   
   #   print(paste(Sys.time(),"no snakes present in",country,":"))
   #   
   # }
   
 }
 
 # write.table(x=CountryExps, file=paste(MyDir, "/Maxent_LUTables/PerCountryExposures.csv",sep=""),
 #             na = "NA", row.names = FALSE, col.names = TRUE, sep=",")
 # 
 
 
 
 
 
 
 ##############################################################
 # show ALL TAXONOMIC GROUPS per COUNTRY that have the highest current exposure index
 # 
 DatSumMod<-read.table(file=paste(MyDir, "/Maxent_LUTables/",projectName,"_SupMatFinFin.csv",sep=""), header=TRUE, sep=",") #load data summary for all species
 DatSumMod<-DatSumMod[DatSumMod$Accepted_Species %in% sppA,]
 region.tab<-read.table(file=paste(MyDir, "/Maxent_LUTables/",projectName,"_RegionData.csv",sep=""), header=TRUE, sep=",") #load data summary for all species
 RegionExps<-read.table(file=paste(MyDir, "/Maxent_LUTables/PerCountryExposures.csv",sep=""), header=TRUE, sep=",") #load data summary for all species
 DatSum<-DatSum[DatSum$Accepted_Species %in% sppA,]
 mygroups<-unique(DatSum$Label)
 
 mycountries<-unique(RegionExps$CountryCode)
 mycountrynames<-unique(RegionExps$CountryName)
 
 for (region in mycountries){
   
   n<- which(mycountries==region)
   s<- which(mycountries==region)
   countryname<-mycountrynames[n]
   ThisRegionExps<-RegionExps[RegionExps$CountryCode==region,]
   
   group.exps<-vector()
   region.groups<-vector()
   
   for (group in mygroups){
     
     m<-which(mygroups==group)
     
     myspp<-DatSum$Accepted_Species[DatSum$Label==group] #list all species in this group
     myrows<-which(ThisRegionExps$Accepted_Species %in% myspp)      #get row numbers of region exposures that belong to species from this group
     
     if(length(myrows)>=1){
       
       mysum<-sum(ThisRegionExps$Exposure_Current[myrows])
       group.exps<-c(group.exps,mysum)
       region.groups<-c(region.groups,group)
       
     }
   }
   
   
   DatSumModRegion<-data.frame(region.groups,group.exps,rep(region,length(region.groups)),rep(countryname,length(region.groups)))
   names(DatSumModRegion)<-c("TaxonomicGroup","Exposure_Current","CountryCode","CountryName")
   DatSumModRegion<-DatSumModRegion[order(DatSumModRegion$Exposure_Current,decreasing=TRUE),]
   
   print(paste(Sys.time(),"exposure for",region,":"))
   print(paste(DatSumModRegion))
   
   if(s==1){
     RegionExpsGroups<-DatSumModRegion
   } else{
     RegionExpsGroups<-rbind(RegionExpsGroups,DatSumModRegion)
   }
   
   # group.exps<-DatSumModRegion$Exposure_Current[1:10]
   # tot.exp<-sum(DatSumModRegion$Exposure_Current,na.rm=TRUE)
   # res.exp<-tot.exp-sum(group.exps,na.rm=TRUE)
   # group.exps<-c(group.exps,res.exp)
   # top10<-DatSumModRegion$MergeShort[1:10]
   # top10<-gsub(" ","\n",c(top10,'other'))
   # 
   # group.exps[11]<-ifelse(group.exps[11] > group.exps[1],group.exps[1],group.exps[11])
   # 
   # png(filename=paste(PlotDir,"/",country,"Main10SppGroups.png",sep=""), units="in", width=16, height=16, res=500) #png smaller file size, A3 for better resolution
   # par(mfrow = c(1,1), mar = c(13,6,6,1), oma=c(0,0.2,0.2,0.2), mgp=c(4,0.5,0))
   # # #mar=c(bottom, left, top, right); oma= outer margins; mpg = axis label locations c(3, 1, 0)
   
   group.exps<-DatSumModRegion$Exposure_Current
   #tot.exp<-sum(DatSumModRegion$Exposure_Current,na.rm=TRUE)
   #res.exp<-tot.exp-sum(group.exps,na.rm=TRUE)
   #group.exps<-c(group.exps,res.exp)
   top10<-DatSumModRegion$TaxonomicGroup
   top10<-gsub(" ","\n",c(top10))
   
   #group.exps[11]<-ifelse(group.exps[11] > group.exps[1],group.exps[1],group.exps[11])
   
   png(filename=paste(PlotDir,"/",region,"_RankedTaxa.png",sep=""), units="in", width=25, height=15, res=500) #png smaller file size, A3 for better resolution
   par(mfrow = c(1,1), mar = c(15,6,6,1), oma=c(0,0.2,0.2,0.2), mgp=c(4,0.5,0))
   
   barplot(group.exps/1000000,las=2,
           col=c(heat.colors(length(top10),1,rev=FALSE)),
           border="grey40",
           space=c(0.1),
           ylim=c(0,max(DatSumModRegion$Exposure_Current/1000000)),
           font.axis=1,
           xlab=NULL,ylab="SHOI [millions]",
           names.arg=top10,
           legend.tex=FALSE,
           main=countryname,
           cex.main=2,cex.axis=1.8,cex.lab=2,cex.names=1.7)
   lines (c(-0.2,length(top10)+1.1),c(0,0), col="black", lwd= 2, lty=1)
   
   # if(group.exps[11] >= group.exps[1]){
   #   
   #   if(round(res.exp/1000000)>=10){
   #     text(x=11.2,y=(group.exps[1]/20/1000000)*16, labels=paste(round(res.exp/1000000)),pos=4, col="black", cex=2)
   #   }else{
   #     text(x=11.3,y=(group.exps[1]/20/1000000)*16, labels=paste(round(res.exp/1000000)),pos=4, col="black", cex=2)
   #   }
   #   
   #   lines (c(11,12.2),c((group.exps[1]/20/1000000)*10,
   #                       (group.exps[1]/20/1000000)*12),
   #          col="black", lwd= 2, lty=1)
   #   lines (c(11,12.2),c((group.exps[1]/20/1000000)*12,
   #                       (group.exps[1]/20/1000000)*14),
   #          col="black", lwd= 2, lty=1)
   #   
   # }
   # 
   
   dev.off()
   
   write.table(x=RegionExpsGroups, file=paste(MyDir, "/Maxent_LUTables/PerCountryExposuresSppGroups_prelim.csv",sep=""),
               na = "NA", row.names = FALSE, col.names = TRUE, sep=",")
   
   # polygon(regionpoly, border = NA,
   #         col=rgb(col2rgb("black")[1],col2rgb("black")[2],col2rgb("black")[3],max=255,alpha=75))
   
 }
 
 
 write.table(x=RegionExpsGroups, file=paste(MyDir, "/Maxent_LUTables/PerCountryExposuresSppGroups.csv",sep=""),
             na = "NA", row.names = FALSE, col.names = TRUE, sep=",")
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

######################################################
# 
# DatSumMod<-read.table(file=paste(MyDir, "/Maxent_LUTables/",projectName,"_SupMatFinFin.csv",sep=""), header=TRUE, sep=",") #load data summary for all species
# 
# DatSumMod$noEDRorSDM<-NA
# DatSumMod$EDRonly<-NA
# DatSumMod$SDMonly<-NA
# DatSumMod$EDRandSDM<-NA
# DatSumMod$SDMsuitOUT<-NA
# DatSumMod$SDMsuitIN<-NA
# DatSumMod$SDMsuitOUTmax<-NA
# 
# #for each species {
# for (sp in spp[which(spp=="Bothrops_alternatus"):length(spp)]){
#   
#   print(paste(Sys.time(),"calculating EDR & SDM overlap for",sp))
#   n<-which(DatSumMod$Accepted_Species==gsub("_"," ",sp))
#   
#   if(n>1){
#     DatSumMod<-read.table(file=paste(MyDir, "/Maxent_LUTables/",projectName,"_SupMatFinFin.csv",sep=""), header=TRUE, sep=",") #load data summary for all species
#   }
#   
#   #read in SDM
#   outfile<-paste(IndSppOutDir,"/",sp,"_Current_CropThreshCDCut.asc",sep="")
#   
#   # unzip raster
#     if(file.exists(paste(outfile,".gz", sep="")) & !file.exists(outfile)){
#     gunzip(filename=paste(outfile,".gz", sep=""),
#            destname=paste(outfile),
#            remove=FALSE, overwrite=TRUE)
#     }
# 
#     # read raster
#     SDM.ras<-raster(paste(outfile), crs="+init=epsg:4326")
# 
#     #read in polygon
#     mypolies <- subset(WHOpolies, WHOpolies$Scientific == gsub("_"," ",sp))
#     #rasterize polygon (10)
#     print(paste(Sys.time(),"cropping snapraster"))
#     snapraster<-crop(GLOBrasterC, extent(SDM.ras))
#     mypolies<-as_Spatial(mypolies)
#     print(paste(Sys.time(),"rasterizing spp polygons"))
#     EDR.ras  <- rasterize(x=mypolies,y=snapraster,as.numeric(10))
#     #make NA=0
#     EDR.ras[is.na(EDR.ras)] <- 0
# 
#     print(paste(Sys.time(),"combining rasters"))
#     
#     #add SDM to polygon 
#     my.ras<-SDM.ras+EDR.ras
#     
#     #remove unzipped file
#     if(file.exists(paste(outfile,".gz", sep="")) & file.exists(outfile)){
#             file.remove(outfile)
#           }
#           if(!file.exists(paste(outfile,".gz", sep="")) & file.exists(outfile)){
#             gzip(filename=paste(outfile))
#           }
# 
#     print(paste(Sys.time(),"extracting values and calculating results"))
#     
#     # get all values
#     myvals <- values(my.ras)
#     myvals<- myvals[!is.na(myvals)]
#     
#     #get all values & calculate length of each (x=0: both say absent, 0<x<1: not in EDR but in SDM, 
#     #                                           x=10: only in EDR but NOT SDM, 10<x<=11: both in EDR & SDM)
#     myvalsNoNo<- myvals[myvals==0]
#     myvalsEDRyes<- myvals[myvals==10]
#     myvalsSDMyes<- myvals[myvals>0 & myvals<=1]
#     myvalsYesYes<- myvals[myvals>10]
#     
#     DatSumMod$noEDRorSDM[n]<-length(myvalsNoNo)
#     DatSumMod$EDRonly[n]<-length(myvalsEDRyes)
#     DatSumMod$SDMonly[n]<-length(myvalsSDMyes)
#     DatSumMod$EDRandSDM[n]<-length(myvalsYesYes)
# 
#     #calculate sum of all 0<x<1 (total suitability predicted outside EDR)
#     #and max (shows if any high suitability predictions lie outside EDR)
#     #and subtract 10 from all 10<x<=11 & calculate sum (total suitability predicted inside EDR)
#     DatSumMod$SDMsuitOUT[n]<-sum(myvalsSDMyes)
#     DatSumMod$SDMsuitOUTmax[n]<-max(myvalsSDMyes,na.rm=TRUE)
#     myvalsYesYes<-myvalsYesYes-10
#     DatSumMod$SDMsuitIN[n]<-sum(myvalsYesYes)
#     
#     #save summary table
#     write.table(x=DatSumMod, file=paste(MyDir, "/Maxent_LUTables/WHOSnakes_SupMatFinFin.csv",sep=""),
#                 na = "NA", row.names = FALSE, col.names = TRUE, sep=",")
#     
#     removeTmpFiles(h=0.1) #removes all temporary files older than 30min, then read hotspot map back in to continue process
#     
# }
# 
# 
# DatSumMod$SDMsuitOUTmax<-ifelse(DatSumMod$SDMsuitOUTmax==(-Inf),0,DatSumMod$SDMsuitOUTmax)
# DatSumMod$SDMtotal<- DatSumMod$SDMonly+DatSumMod$EDRandSDM
# DatSumMod$SDMtotal<-ifelse(DatSumMod$SDMtotal<=0,0.01,DatSumMod$SDMtotal)
# DatSumMod$SDMtotal<-ifelse(DatSumMod$SDMtotal< DatSumMod$EDRonly,DatSumMod$EDRonly,DatSumMod$SDMtotal)
# 
# 
# 
# DatSumMod$SDMsuittotal<-DatSumMod$SDMsuitOUT+DatSumMod$SDMsuitIN
# DatSumMod$SDMsuittotal<-ifelse(DatSumMod$SDMsuittotal<=0,0.01,DatSumMod$SDMsuittotal)
# 
# DatSumMod$SPrecs<-ifelse(DatSumMod$SPrecs==(-Inf),0,DatSumMod$SPrecs)
# 
# 
# #save summary table
# write.table(x=DatSumMod, file=paste(MyDir, "/Maxent_LUTables/WHOSnakes_SupMatFinFin.csv",sep=""),
#             na = "NA", row.names = FALSE, col.names = TRUE, sep=",")
# #write.table(x=DatSumMod, file=paste(MyDir, "/2_Maxent_LUTables/WHOSnakes_SupMatFinFin.csv",sep=""),
# #            na = "NA", row.names = FALSE, col.names = TRUE, sep=",")




# 
# 
# 
# 
# DatSumMod<-read.table(file=paste(MyDir, "/Maxent_LUTables/",projectName,"_SupMatFinFin.csv",sep=""), header=TRUE, sep=",") #load data summary for all species
# 
# print(paste(Sys.time(),"plotting results across all species"))
# 
# png(filename=paste(PlotDir,"/EDRSDMboxplots.png",sep=""), units="in", width=8, height=16, res=500) #png smaller file size, A3 for better resolution
# par(mfrow = c(2,1), mar = c(10,10,5,1), oma=c(0,0.5,2,0.5), mgp=c(5,0.5,0))
# # #mar=c(bottom, left, top, right); oma= outer margins; mpg = axis label locations c(3, 1, 0)
# 
# #1. plot boxplot of areas (length) inside and outside EDR
# boxplot(list(
#   DatSumMod$SDMonly/DatSumMod$SDMtotal,
#   DatSumMod$EDRandSDM/DatSumMod$SDMtotal,
#   DatSumMod$EDRonly/DatSumMod$SDMtotal,
#   DatSumMod$SDMsuitOUT/DatSumMod$SDMsuittotal,
#   DatSumMod$SDMsuitIN/DatSumMod$SDMsuittotal,
#   DatSumMod$SDMsuitOUT/DatSumMod$SDMonly,
#   DatSumMod$SDMsuitIN/DatSumMod$EDRandSDM),
#   ylab=paste("Fraction of total SDM value (0-1)"),
#   names=c("Fraction\nCommission\nArea","Fraction\nCongruence\nArea","Fraction\nOmission\nArea", 
#           "Fraction\nCommission\nSuitability","Fraction\nCongruence\nSuitability",
#           "Average\nCommission\nSuitability","Average\nCongruence\nSuitability"),
#   col=c("white","grey90","grey50","grey30","grey10")
# )
# 
# #commission=false positive, omission=false negative
# dev.off()
# 
# 
# # png(filename=paste(PlotDir,"/EDRSDMscatterplots.png",sep=""), units="in", width=10, height=20, res=500) #png smaller file size, A3 for better resolution
# # par(mfrow = c(3,1), mar = c(5,5,1,1), oma=c(1,0.5,1,0.5), mgp=c(2,0.5,0), bty="l")
# # #mar=c(bottom, left, top, right); oma= outer margins; mpg = axis label locations c(3, 1, 0)
# 
# # plot(log(DatSumMod$SPrecs+1),DatSumMod$SDMtotal/DatSumMod$SDMtotal,col="white",ylim=c(0,1))
# # points(log(DatSumMod$SPrecs+1),DatSumMod$SDMonly/DatSumMod$SDMtotal,col="red")
# # points(log(DatSumMod$SPrecs+1),DatSumMod$EDRonly/DatSumMod$SDMtotal,col="grey")
# # points(log(DatSumMod$SPrecs+1),DatSumMod$EDRandSDM/DatSumMod$SDMtotal,col="orange")
# # abline(lm((DatSumMod$SDMonly/DatSumMod$SDMtotal) ~ log(DatSumMod$SPrecs+1)), col="red")
# # abline(lm((DatSumMod$EDRonly/DatSumMod$SDMtotal) ~ log(DatSumMod$SPrecs+1)), col="grey")
# # abline(lm((DatSumMod$EDRandSDM/DatSumMod$SDMtotal) ~ log(DatSumMod$SPrecs+1)), col="orange")
# # 
# # plot(log(DatSumMod$EDRonly),DatSumMod$SDMtotal/DatSumMod$SDMtotal,col="black",ylim=c(0,1))
# # points(log(DatSumMod$EDRonly),DatSumMod$SDMonly/DatSumMod$SDMtotal,col="red")
# # #points(DatSumMod$EDRonly,DatSumMod$EDRonly)
# # points(log(DatSumMod$EDRonly),DatSumMod$EDRandSDM/DatSumMod$SDMtotal,col="orange")
# 
# # plot(log(DatSumMod$SPrecs+1),DatSumMod$SDMsuittotal/DatSumMod$SDMsuittotal,pch=16,col="white",ylim=c(0,1))
# # points(log(DatSumMod$SPrecs+1),DatSumMod$SDMsuitOUT/DatSumMod$SDMsuittotal,col="grey",pch=16,cex=0.7)
# # points(log(DatSumMod$SPrecs+1),DatSumMod$SDMsuitIN/DatSumMod$SDMsuittotal,col="black",pch=8,cex=0.7)
# # abline(lm((DatSumMod$SDMsuitOUT/DatSumMod$SDMsuittotal) ~ log(DatSumMod$SPrecs+1)), col="grey")
# # abline(lm((DatSumMod$SDMsuitIN/DatSumMod$SDMsuittotal) ~ log(DatSumMod$SPrecs+1)), col="black")
# # cor(x=log(DatSumMod$SPrecs+1),y=(DatSumMod$SDMsuitOUT/DatSumMod$SDMsuittotal),method = c("spearman"))
# # cor(x=log(DatSumMod$SPrecs+1),y=(DatSumMod$SDMsuitIN/DatSumMod$SDMsuittota),method = c("spearman"))
# 
# png(filename=paste(PlotDir,"/EDRSDMcongruence.png",sep=""), units="in", width=10, height=20, res=500) #png smaller file size, A3 for better resolution
# par(mfrow = c(3,1), mar = c(7,7,1,1), oma=c(1,1,1,0.5), mgp=c(5,2,0), bty="l")
# 
# boxplot((DatSumMod$SDMsuitOUT/DatSumMod$SDMsuittotal)*100~round(log(DatSumMod$SPrecs+1),0),
#         xlim=c(1.8,10.2),col="grey",varwidth=TRUE,pch=8,box=FALSE,notch=FALSE,
#         ylab="Percent suitability commission",xlab=NULL,cex.lab=2,cex.axis=2)
# 
# 
# boxplot((DatSumMod$SDMsuitIN/DatSumMod$SDMsuittotal)*100~round(log(DatSumMod$SPrecs+1),0),
#         xlim=c(1.8,10.2),col="grey",varwidth=TRUE,pch=8,box=FALSE,notch=FALSE,
#         ylab="Percent suitability congruence",xlab=NULL,cex.lab=2,cex.axis=2)
# 
# boxplot((DatSumMod$EDRonly/(DatSumMod$EDRandSDM+DatSumMod$EDRonly))*100~round(log(DatSumMod$SPrecs+1),0),
#         col="grey",pch=8,lty=1,varwidth=TRUE,xlim=c(1.8,10.2),
#         ylab="Percent distribution omission",xlab="Log number of occurrence records",cex.lab=2,cex.axis=2)
# 
# dev.off()
# 
# 
# 
# 
# # plot(log(DatSumMod$EDRonly),DatSumMod$SDMsuittotal/DatSumMod$SDMsuittotal,col="black",ylim=c(0,1))
# # points(log(DatSumMod$EDRonly),DatSumMod$SDMsuitOUT/DatSumMod$SDMsuittotal,col="red")
# # points(log(DatSumMod$EDRonly),DatSumMod$SDMsuitIN/DatSumMod$SDMsuittotal,col="orange")
# 
# 
# 
# 
# #############################################################################################
# #############################################################################################
# #############################################################################################